Conditional Mean and Var of GARCH

• Conditional Mean: \(0\)

• Conditional Variance : \(\sigma_t^2\)

Ex: Daily SPY

Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014

  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/SPY.csv", header=T)

## Warning in if (!header) rlabp <- FALSE: the condition has length > 1 and
## only the first element will be used

## Warning in if (header) {: the condition has length > 1 and only the first
## element will be used

  source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
    
  head(D)

##   ixDay      Date   Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1     1  1/3/2000    Monday 1   148.25   148.25  143.88    145.44
## 2     2  1/4/2000   Tuesday 2   143.53   144.06  139.64    139.75
## 3     3  1/5/2000 Wednesday 3   139.94   141.53  137.25    140.00
## 4     4  1/6/2000  Thursday 4   139.62   141.50  137.75    137.75
## 5     5  1/7/2000    Friday 5   140.31   145.75  140.06    145.75
## 6     6 1/10/2000    Monday 6   146.25   146.91  145.03    146.25
##   SPY.Volume SPY.Adjusted
## 1    8164300       110.33
## 2    8089800       106.01
## 3   12177900       106.20
## 4    6227200       104.50
## 5    8066500       110.57
## 6    5741700       110.94

  X <- ts(D[,7], start=1)
  plot(X)

  X1 <- diff(log(X))
  plot(X1)

  layout(matrix(1:2, 1, 2))
  acf(X1)
  acf(X1^2)

  layout(1)

  # Fit GARCH model 
  library(fGarch)
  Fit1 <-  garchFit(~ garch(1,1), data=X1, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    
  ##  Fit1@fit$par          #-- estimated parameters
  ##  Fit1@residuals        #-- this is not the garch residuals!!!! same as X1
  ##  Fit1@sigma.t          #-- estimated sig_t
  ##  Fit1@residuals/Fit1@sigma.t       #-- this is the (standardized) GARCH residuals

  print(Fit1@fit$ics)   #-- AIC and BIC are here

##       AIC       BIC       SIC      HQIC 
## -6.470319 -6.465360 -6.470321 -6.468556

  plot(X1)
  lines( 1.96*Fit1@sigma.t, type="l", col="red")
  lines(-1.96*Fit1@sigma.t, type="l", col="red")

  res1  <- X1/Fit1@sigma.t   #-- this is the residuals
  Randomness.tests(res1)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##      BL15 BL20 BL25  ML15  ML20 JB    SD
## [1,]    0    0    0 0.719 0.855  0 1.001

Write \(\sigma^2_t\) using \(Y_t^2\)

• Starting with the definition of GARCH(1,1),

  \[ \begin{align} \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \\\\ &= \omega + \alpha Y_{t-1}^2 + \beta \Big(\omega + \alpha Y_{t-2}^2 + \beta \sigma^2_{t-2} \Big)\\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \sigma^2_{t-2} \\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \Big(\omega + \alpha Y_{t-3}^2 + \beta \sigma^2_{t-3} \Big) \\ \\ &= \omega + \beta \omega + \beta^2 \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \alpha Y_{t-3}^2 + \beta^3 \sigma^2_{t-3} \end{align} \]

• Continuing, we get \[ \begin{align} &= \omega (1 + \beta + \beta^2 + \cdots ) + \alpha \sum_{i=0}^k \beta^i Y_{t-1-i}^2 + \beta^{k+1} \sigma^2_{t-1-k} \\\\ &= \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \end{align} \]

• That means if we keep going, we can write \[ \sigma_t^2 \hspace3mm = \hspace3mm \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \]

• We will use the truncated and estimated version of this, \[ \hat \sigma_t^2 \hspace3mm = \hspace3mm \frac{\hat \omega}{1- \hat \beta} \hspace3mm + \hspace3mm \hat \alpha \sum_{i=0}^{t-1} \hat \beta^i Y_{t-1-i}^2 \]

Quasi-MLE

• Quasi-Maximum Likelihood Estimator, means we will use MLE for normal errors, even in the case of non-normality. Theoretically, for GARCH parameters, it is shown that Normal MLE is still consistent under non-normal errors.

  x <- seq(-10,10, .1)
    
  #- Standardized t-distribution 
  plot(x, dstd(x, mean = 0, sd = 1, nu = 5), type="l" )
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )  

  #- Generalized Error Distribution
  plot(x, dged(x, mean = 0, sd = 1, nu = 4), type="l" )
  lines(x, dnorm(x, mean = 0, sd=1) , col="red" )

  #- Skewed-Standardized t-distribution 
  plot( x, dsstd(x, mean = 0, sd = 1, nu = 5, xi = 1.5), type="l" )
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

  #- Skewed-Generalized Error Distribution
  plot(x, dsged(x, mean = 0, sd = 1, nu = 5, xi = 1.5) , type="l" )
  lines(x, dnorm(x, mean = 0, sd = 1), col="red" )

Simulation Study

True parameter (.024, .1, .8)

• True conditional distribution: Normal Estimated using: Normal

• True conditional distribution: std(5) Estimated using: Normal

• True conditional distribution: sged(skew=.7, shape=1.45) Estimated using: Normal

  #--- Simulation showing how error distribution affect parameter estimation in GARCH(1,1) ---
  n   <- 1000
  itt <- 1000
  theta <- c(.024, .1, .8)
  
  Param11 <- matrix(rep(0,3*itt), itt, 3)
  for (i in 1:itt) {
  
    spec <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="norm")
    x    <- garchSim(spec, n = n, extended=FALSE)
  
    est2 <-  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    Param11[i,] <- coef(est2)
  }
  
  Param12 <- matrix(rep(0,3*itt), itt, 3)
  for (i in 1:itt) {
  
    spec <- garchSpec(model = list(shape=5, omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="std")
    x     <- garchSim(spec, n = n, extended=FALSE)
  
    est2 <-  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    Param12[i,] <- coef(est2)
  }

  Param13 <- matrix(rep(0,3*itt), itt, 3)
  for (i in 1:itt) {
  
    spec <- garchSpec(model = list(skew=.9,  shape=1.45, omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="sged")
    x    <- garchSim(spec, n = n, extended=FALSE)
  
    est2 <-  garchFit(~ garch(1,1), data=x, cond.dist="norm",  include.mean = FALSE, trace = FALSE)
    Param13[i,] <- coef(est2)  
  }

  layout(matrix(1:9, 3, 3, byrow=T))
  hist(Param11[,1], 10, xlim=c(0,.5)); hist(Param11[,2], 10, xlim=c(0,.5)); hist(Param11[,3], xlim=c(0.4,1), breaks=10)
  hist(Param12[,1], 10, xlim=c(0,.5)); hist(Param12[,2], 10, xlim=c(0,.5)); hist(Param12[,3], xlim=c(0.4,1), breaks=10)
  hist(Param13[,1], 10, xlim=c(0,.5)); hist(Param13[,2], 10, xlim=c(0,.5)); hist(Param13[,3], xlim=c(0.4,1), breaks=10)

  P11 <- (Param11 - t(matrix(theta, 3, itt) ))^2
  P12 <- (Param12 - t(matrix(theta, 3, itt) ))^2
  P13 <- (Param13 - t(matrix(theta, 3, itt) ))^2
  
  apply(P11, 2, mean)
  apply(P12, 2, mean)
  apply(P13, 2, mean)
  
  apply(P11, 2, mean)/apply(P11, 2, mean)
  apply(P12, 2, mean)/apply(P11, 2, mean)
  apply(P13, 2, mean)/apply(P11, 2, mean)

MSE1/MSE1   #- Relative MSE 
[1] 1 1 1

MSE2/MSE1
[1] 1.584734 3.396001 1.831015

MSE3/MSE1
[1] 1.198293 1.613116 1.355613

Stationarity Cond’n for GARCH

• GARCH(1,1) is weakly stationary if \[ E( \log(\beta + \alpha e_t^2)) < 0 \]

• This is satisfied if \[ \alpha + \beta < 1. \]

Using Simulated Data

• See how closely true sigma is estimated and predicted, using simulated GARCH

library(quantmod)
library(fGarch)
source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt") 

  theta = c(.5, .7, .2)
  spec1 <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]),  cond.dist="norm")
    
  #- Generate GARCH
  set.seed(66343)
  X <- garchSim(spec1, n = 1000, extended=T)   #extended=T:  now  Y has 3 columns: garch sigma eps

## Warning in if (extended) ans else ans[, "garch"]: the condition has length
## > 1 and only the first element will be used

  X <- as.xts(X)
  X <- xts(X, order.by=as.Date(index(X)))

  Y        <- X[,1]    #- simulated GARCH
  sig.true <- X[,2]    #- this is true sigma

  plot(Y, main="Simulated GARCH")

  plot(sig.true, main="True Sigma_t")

  #- Estimate GARCH
  out1  <-  garchFit(~ garch(1,1), data=Y, cond.dist="norm", include.mean = FALSE, trace = FALSE)
  sig.estim <- xts(out1@sigma.t, order.by=index(Y))


  #- compare true sigma vs estimated sigma
  head( cbind(sig.true, sig.estim) )

##                sigma       ..2
## 2015-07-14 0.9844680 1.7773189
## 2015-07-15 0.8719402 1.0198346
## 2015-07-16 0.8528739 0.9322072
## 2015-07-17 0.9465185 1.0167695
## 2015-07-18 0.9430606 1.0033031
## 2015-07-19 1.1075998 1.1624516

  plot(cbind(sig.estim, sig.true), col=c("red", "black"), main="True Sigma (Blk) vs Estimated (Red)")

  #- Compute the sigma for tomorrow (forecast sigma_t)
  w <- out1@fit$par[1] 
  a <- out1@fit$par[2] 
  b <- out1@fit$par[3] 
  sig.pred <- xts( sqrt( w + a*Y[1000]^2 + b*out1@sigma.t[1000]^2 ), order.by=index(Y[1000])+1)

  #- This does the same thing (10-day forecast sigma_t)
  sig.pred2 <- predict(out1)   
  sig.pred2 

##    meanForecast meanError standardDeviation
## 1             0  1.061806          1.061806
## 2             0  1.267963          1.267963
## 3             0  1.411732          1.411732
## 4             0  1.517237          1.517237
## 5             0  1.596835          1.596835
## 6             0  1.657922          1.657922
## 7             0  1.705340          1.705340
## 8             0  1.742442          1.742442
## 9             0  1.771642          1.771642
## 10            0  1.794720          1.794720

  sig.pred2[1,3] 

## [1] 1.061806

  #- Rolling 1-day prediction of sig.t
  sig.pred1 <- numeric(0)
  date.stp <- numeric(0)
  for (i in 1:750){
    T <- Y[(1:250)+(i-1)]
    n <- length(T)
    out1  <-  garchFit(~ garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
    sig.pred1[i] <- predict(out1)[1,3]   #- sig.t prediction
    date.stp[i]  <- index(T[n])+1
  }
  sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))


  #- Plot and compare estimated vs rolling predicted
  plot(cbind(sig.pred, sig.estim, sig.true),
    col=c("blue", "red", "black"),
    main="sigma.t:  True (Blk) vs Estim (Red) vs Pred(Blue)")

  plot(cbind(sig.pred, sig.estim, sig.true)['2016'],
    col=c("blue", "red", "black"),
    main="sigma.t:  True (Blk) vs Estim (Red) vs Pred(Blue)")

  plot(cbind(sig.pred, sig.estim, sig.true)['2016-04'],
    col=c("blue", "red", "black"),
    main="sigma.t:  True (Blk) vs Estim (Red) vs Pred(Blue)")

  ts.plot(cbind(Y, 1.96*sig.true, -1.96*sig.true, 1.96*sig.pred, -1.96*sig.pred),
      col=c("black", "black", "black", "blue", "blue"))

  ts.plot(cbind(Y, 1.96*sig.true, -1.96*sig.true, 1.96*sig.pred, -1.96*sig.pred)['2016-06::2016-12'],
      col=c("black", "black", "black", "blue", "blue"))

Using SPX

• See how close estimated sigma and rolling predicted sigma, using real data.

  getSymbols("^GSPC")     #- download S\&P500 index from Yahoo

## [1] "GSPC"

  SPX <- Ad(GSPC)
  plot(SPX)

  plot(diff(log(SPX)))

  #- we model this with GARCH (this time we don't know the true sigma.t)
  Y <- diff(log(SPX["2013::2016"]))[-1]
  length(Y)

## [1] 1007

  #- Estimate GARCH
  out1  <-  garchFit(~ garch(1,1), data=Y, cond.dist="norm", include.mean = FALSE, trace = FALSE)
  sig.estim <- xts(out1@sigma.t, order.by=index(Y))

  #- Rolling 1-day prediction of sig.t
  sig.pred1 <- numeric(0)
  date.stp  <- numeric(0)
  for (i in 1:750){
    T <- Y[(1:250)+(i-1)]
    out1  <-  garchFit(~ garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
    sig.pred1[i] <- predict(out1)[1,3]   #- sig.t prediction
    date.stp[i]  <- index(T[length(T)]+1)
  }
  sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))

  #- Compare estimated sigma vs predicted sigma - 
  plot(cbind(sig.pred, sig.estim), col=c("blue", "red"))

  #- Compare daily CI from two sigmas
  plot(as.numeric(Y["2016"]), type="h")
  lines(as.numeric( 1.65*sig.estim["2016"]), col="red")
  lines(as.numeric(-1.65*sig.estim["2016"]), col="red")
  lines(as.numeric( 1.65*sig.pred["2016"]),  col="blue")
  lines(as.numeric(-1.65*sig.pred["2016"]),  col="blue")

More Examples

• Cryer : CREF stock values

• Shumway : NYSE us GNP

• Cowpertwait : SP500, Southern hem Temp.